Dynamics of quantum coherences at strong coupling to a heat bath 
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The standard approach for path integral Monte Carlo simulations of open quantum systems is 
extended as an efficient tool to monitor the time evolution of coherences (off-diagonal elements of the 
reduced density matrix) also for strong coupling to environments. Specific simulations are performed 
for two level systems embedded in Ohmic and sub-Ohmic reservoirs in the domains of coherent and 
incoherent dynamics of the polarization. In the latter regime, the notorious difficulty to access the 
long time regime is overcome by combining simulations on moderate time scales with iteratively 
calculated initial densities. This allows to extract relaxation rates for sub-Ohmic environments 
at finite temperatures and over a broad range of couplings and to compare them to analytical 
predictions. The time evolution of the von Neumann entropy provides insight into the quantum 
phase transition at thermal equilibrium from a delocalized to a localized state at zero temperature. 

PACS numbers: 03.65.Yz, 05.10.Ln, 03.67.-a, 73.63.-b 
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I. INTRODUCTION 

The dynamics of open quantum systems has again re- 
ceived considerable interest in the last years 1 ^, mainly 
due to the progress in fabricating and manipulating quan- 
tum devices with unprecedented accuracy. Examples in- 
clude complex superconducting circuits and tailored set- 
ups for ensembles of ultra-cold atomic gases. In fact, 
the interaction to reservoirs has even been found as a re- 
source to drive the time evolution of complex system into 
desired final states^. 

Accordingly, there have been substantial efforts to de- 
velop new and to extend existing numerical method- 
ologies to treat open system dynamics beyond conven- 
tional perturbative approaches like e.g. master equa- 
tions. Among them are those that work with the re- 
duced density operator of the relevant system and those 
that treat the full Hilbert space of system + reservoir. 
The latter comprise renormalization group methods 4 
(NRG/DMRG) and advanced basis set techniques 5 
(MCTDH), the former ones include the quasi-adiabatic 
propagator approach^ (QUAPI), stochastic Schrodinger 7 ' 
and Liouville-von Neumann 8 (SLN) formulations, and 
path integral Monte Carlo schmemes (PIMC)£ii£. Each 
of these methods has its strengths and its limitations. 
The most challenging situations appear in the domains 
of strong coupling to heat baths with non-Ohmic spectral 
densities at low but elevated temperatures. There, the 
non-Markovian retardation is influential and system-bath 
correlations are affected by both quantum and thermal 
fluctuations. 

The PIMC approach belongs to the very few methods 
that yield numerically exact results for the dynamics in 
all ranges of parameter space, particularly, in the deep 
non-perturbative regime. It is based on the formally ex- 
act Feynman- Vernon formulation of the reduced dynam- 
ics in terms of path integrals^. The impact of the reser- 
voir is captured by a so-called influence functional which 
contains the force-force correlation of the bath and thus 
carries a retarded interaction in time between system 



paths. Impressive progress has been achieved in recent 
years to extend it to a very efficient and powerful tool for 
the quantum dynamics of systems with discrete Hilbert 
space (tight-binding systems)^. Single and correlated 
charge transfer along molecular chain o 10 ' 11 , the impact 
of external driving fields 1 ^ and steady states in charge 
transfer through quantum dots and molecular contacts 1 ^ 
have been studied. 

While in these applications the reservoirs are basi- 
cally Ohmic (typically with a large cut-off), recently sub- 
Ohmic reservoirs and their entanglement with the system 
of interest gained attention 1 ^—. These reservoirs have 
spectral densities of the form J(w) ~ au s with typi- 
cal coupling constant a and spectral exponent s. They 
are not only generic models to better understand fun- 
damental concepts in quantum mechanics, but also to 
explain at least qualitatively experimental observations 
for superconducting qubits 1 ^ and quantum dots^, quan- 
tum impurity systems^ 1 -, nanomechanical oscillators^, 
and trapped ion systems^. In this context, a generic 
model is that of a two level system (artificial spin- i) in- 
teracting with a heat bath (spin-boson model 1,14 , SBM). 
The SBM has been analyzed in variety of arrangements in 
the past 1 -, however, its sub-Ohmic properties have been 
left basically untouched mainly due to the lack of pow- 
erful numerical schemes. In the last years, this gap has 
been closed by applying advanced techniques mentioned 
above. It turns out that the dominance of low frequency 
modes causes strong retardation effects in the dynam- 
ics, quantum phase transitions in thermal equilibrium at 
zero temperature 1 - 5 - accompanied by substantial system- 
reservoir entanglement^. Very recently, the domain of 
strong friction, unaccessible by alternative approaches, 
has been investigated with PIMC simulations and re- 
vealed the persistence of coherences^ 5 -. 

For that purpose, the approach has to be extended to 
measure not only polarizations (populations), but also 
off-diagonal elements of the density. In this work we go 
one step further and use the knowledge of the full reduced 
density matrix to improve the PIMC technique to cover 
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also long times. Namely, due to the dynamical sign prob- 
lem, i.e. an exponentially decreasing signal to noise ratio 
with growing simulation time, for the PIMC method are 
not reachable using the PIMC method. The idea is now 
to glue together sequences of PIMC simulations for mod- 
erate times with successively determined initial densities. 
In general, this procedure does not work due to intricate 
system-bath correlations. For moderate to strong friction 
though, coherences in the system dynamics are strongly 
suppressed even at very low temperatures and the re- 
laxation towards thermal equilibrium is monotone. Cor- 
responding rates are still quantum mechanical and can 
be extracted efficiently from relatively cheap long time 
simulations based on this concept. 

The paper is organized as follows: The model and rel- 
evant observables are introduced in Sec. [TT] to lay the 
basis for a brief description of the PIMC in Sec. Mil and 
analytical treatments in Sec. IIVI The implementation of 
the measurement operator for the coherences and the ex- 
tended PIMC is first studied in Sec. [V] for the much sim- 
pler Ohmic case as a test-bed. In Sec. lVII the challenging 
sub-Ohmic SBM is considered, particularly its relaxation 
rates in the incoherent domain and the entropy. Conclu- 
sions are given in the end. 



II. MODEL AND OBSERVABLES 

The spin-boson model (SBM) describes a two-level sys- 
tem (TLS) linearly coupled to a heat bath environment. 
This is modeled as an ensemble of harmonic oscillators 
according to the rationale that reservoirs with macro- 
scopic many degrees of freedom carry Gaussian fluctua- 
tion properties. Here, we focus on the symmetric situa- 
tion 
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with degenerate eigenstates ±1) of the Pauli a z operator 
coupled by a tunneling energy h/S.. In the continuum 
limit, the distribution of the bath frequencies is described 
by a spectral density J(lu) — ^J2 a ^{^a — w). 

We are interested in the time evolution of the observ- 
ables 

P„(t) = (a v (t)) 

= Tr s {a v p(t)} , u = x,y,z (2) 

which fully characterize the reduced density 

p(t) = Tr B {exp(-iHt/h)W(0) ^{iHt/h)} . (3) 

The initial state is assumed to be a product state of the 
form 
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with the collective bath operator E = ^ Q hXaip^ + b a ), 
the reservoir's partition function Zb, and the initial state 
po of the bare TLS. The polarization parameter p allows 
to start with shifted initial bath distributions. For exam- 
ple, for p = ±1 the reservoir is equilibrated to an initial 
spin state | ± 1) (polarized), while for p = one regains 
the bare heat bath. 

In contrast to previous PIMC simulations, we consider 
for the initial state po not only eigenstates of the a z - 
operator (localized spin states), but also generalized ini- 
tial preparations 



V v=x,y,z J 



(5) 



with properly chosen initial values Pv(0) fulfilling the 
constraint 'Yj V P v {Q) 2 = 1- 



III. PATH INTEGRAL MONTE CARLO 
TECHNIQUE 

A non-perturbative treatment of the dynamics ^ is 
obtained within the path integral formulation. The P v is 
expressed along a Keldysh contour with forward a and 
backward a' paths 1 . The impact of the environment ap- 
pears as an influence functional introducing arbitrarily 
long-ranged interactions between the paths. Switching 
to sum and difference paths, respectively, 
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one arrives at 
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with the contribution A v in absence of dissipation includ- 
ing the measurement operator and the influence func- 
tional 

*[C,v] = fdt' f dt"0')\Q'(t'-t"X(t") 
Jo Jo 1 

+ iQ"(t' -t")r)(t")] -ip [ dt'((t')Q'(t'), 
J Jo 

(8) 

where the last term accounts for the initial preparation of 
the bath2&. Here, one sums over all paths which connect 
77(0), C(0) with r](t),((t) where the initial values have to 
be chosen according to the initial preparation ([5]) and the 
final endpoints according to the measured observable © . 
The kernel Q = Q'+iQ" is related to the bath correlation 
Q(t) = (£(t)£(0))/h 2 and reads 



Q(t) 
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( 1 — cos ut ) + i sin cot 
(9) 
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The most efficient methodology to evaluate the above 
path integrals non-perturbatively is the path integral 
Monte Carlo approach (PIMC)*k^. The procedure starts 
from a Trotter-Suzuki discretization of the total time in- 
terval t with time increment r = t/q and Trotter number 
q. The full expression ([7]) reads with PI = P u ((i — l)r) 



H = £ Po(r? l5 Ci) n^(Ci,o+i) 



3 = 1 

xM Uti [C,rj\ exp(-$[C,77])- 
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(10) 



with matrices KK, ni _ ni+1 containing the free propagation 
of the TLS and M v ^ being the measurement operator for 
the observable a v at the time step [i — l)r (for details 
seei^ and the Appendix). 

The impact of the reservoir is completely determined 
by the spectral density. A common form is given by 



J s (w) = 2-Kauj], s uj s exp(— ui/oj c ) 



(11) 



where a is a dimensionless coupling strength, uj c a soft 
ultra- vilolett-cutoff, and s the spectral exponent. For 
s = 1, the bath is called Ohmic, while for < s < 1 
it is referred to as sub-Ohmic. In both cases, analytical 
expressions for Q(t) are known. 

The numerical evaluation of the expression (fTU|) by 
means of Monte Carlo methods is hampered by the so- 
called dynamical sign problem which originates from the 
interference structure of quantum mechanics and leads 
to an exponential decrease of the signal to noise ratio for 
longer times. In the past, several techniques have been 
developed to substantially soothe its impact on PIMC 
simulation*^. The general idea is to reduce the sam- 
pling space by integrating out exactly large parts of the 
configuration space, i.e., to sum over the paths r\ in (|10[) . 
Only the remaining q-fold sum over the relative coor- 
dinates Ci<i<<? i s performed by a Metropolis sampling. 
This way, simulations yield numerically exact results up 
to times where equilibration sets in with error bars of typ- 
ically less than 1%. We will discuss below a procedure 
how to further increase the simulation time by linking 
segments of PIMC on moderate time scales. 



IV. ANALYTICAL TREATMENT 

An approximate and well-studied treatment of the 
SBM which captures also strong coupling and low 
temperatures is the noninteracting blip approximation 
(NIBA) 1 . It is based on a separation of scales be- 
tween times where p(t) resides in a diagonal state (so- 
journ) and those where it is in an off-diagonal state 
(blip). In the NIBA only reservoir induced intra-blip 
and blip-proceeding sojourn interactions in the influence 
functional (jHJl are kept, while remaining correlations are 
dropped. Although this approximation is expected to fail 
in regions of parameter space with long memory effects, 



it provides in many cases quantitatively accurate results 
and in others at least a qualitative description of the 
dynamics. For further details we refer to the extensive 
literatures^. 

For the time dependent polarization P z (t) — (a z (t)) 
with the initial preparation P z (0) = fi = 1 one derives 



dP z (t) 



dt 



du K(t - u)P z {u) 



with kernel 



K(t) = A 2 e- Q ' (t) cos[Q"(t)}. 



(12) 
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In ranges of parameters space where K(t) is sufficiently 
short ranged, the dynamics (|12l) can be approximated by 



dP z {t) 
dt 



-Fniba Pz{t) 



(14) 



with a rate Tniba = J °° du K(u) describing monotone 
decay. 

Further, for the imaginary part of the coherences one 
has P y {t) = —P z (t)/A, while the real part P x follows for 
a symmetric system from 



P x (t) = A [ du e' Q ' {u} sinQ'» . 
Jo 
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As mentioned above, a major goal of this work is to 
explore the options arising when there is also access to 
off-diagonal elements of the reduced density operator and 
when using generalized initial states ([5]). We start with 
some results for an Ohmic environment as a test-bed for 
the more challenging sub-Ohmic case considered after- 
wards. First, PIMC data for P x are presented and com- 
pared to NIBA predictions. This provides the basis for a 
new technique to extend the simulation time for PIMC. 



A. Off-diagonal elements 

It is well-known that the Ohmic SBM displays damped 
oscillatory dynamics in P z (t) for weaker couplings and/or 
lower temperatures (coherent domain) and incoherent de- 
cay for stronger friction and/or higher temperat ure a 1 1 14 . 
At zero temperature, this transition occurs at a = 1/2 
for lo c ^> A. These and other properties have been stud- 
ied numerically by several approaches including PIMC 
and also in comparison with NIBA. Less attention has 
been paid to P x though, for which it is believed that 
NIBA predictions are less reliable. For example, in the 
incoherent regime, Eq. (fT5j) predicts at T = a finite 
value 
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(16) 
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At 

FIG. 1: Dynamics of P x (t) for Ohmic spectral density at T — 
0. PIMC results (dots with error bars) and analytical NIBA 
results (black lines) for various values of the coupling a and 
the cut-off cj c . 



for a > 1/2 while it diverges for a < 1/2. 

The implementation of a x as measurement operator in 
the PIMC is not trivial. Namely, the approach is only 
efficient when one extracts observables at intermediate 
times u from a single run over the total time t > u. 
Hence, since the influence functional is non-local in time, 
the effective measurement operator in the PIMC corre- 
sponding to o~ x must also include contributions from the 
bath 25 (see also Appendix). 

Fig. [T]shows a comparison between PIMC data and the 
NIBA at zero temperature. A polarized initial prepara- 
tion is chosen with P z (0) = 1. For longer times, P x (t) in- 
deed increases stronger for smaller values of the coupling, 
while it quickly saturates for stronger damping. In this 
latter range (a > 0.8), NIBA is in excellent agreement 
with the data. For 0.8 > a > 0.5, i.e., in the incoherent 
domain, it still provides qualitatively correct results over 
the time period considered, but fails completely in the 
coherent range a < 0.5. Further, the larger the cut-off, 
the better is the agreement with the PIMC simulations. 
Increasing a and uj c tends to suppress coherences and 
thus pushes the system towards a localized state which 
exists for a > 1. 

For finite temperatures T > 0, the dynamics of P x (t) 
at moderate damping remains qualitatively the same (not 
shown) and is strongly suppressed at high temperatures 
when Q' dominates against Q" in (|T5|) . 



B. Long-time simulations: CH-PIMC 

The time dependent density operator of the full com- 
pound obeys the semi-group property 

W(t) = U{t,t s )W{t s )U\t,t s ) 
= U(t,0)W(0)U f (t,0) 

(17) 



with U(t,t') — exp[—iH(t — t')/K). In general, this is 
not true for the reduced density operator ©. How- 
ever, there are two domains in parameter space where 
this symmetry holds at least approximately. One is the 
domain of very weak coupling, sufficiently high temper- 
atures, and large cut-off frequencies, i.e., aAhfi <C 1 
and A <C w c . Then, a perturbative treatment including 
a coarse graining in time (Born-Markov approximation) 
leads to master equations for p(t) which are local in time 
and are determined by a time-independent generator. In 
the other domain, coherences (off-diagonal elements) in 
p(t) are suppressed due to the bath such that effectively 
the full density operator at intermediate times takes the 
form W(t s ) » p(t s ) exp[-/3(iT B - p s £)]/Z B . This can 
be interpreted as a quasi-classical regime, where, how- 
ever, temperature does not need to be high, but rather 
friction sufficiently strong. Accordingly, quantum effects 
may still substantially influence the relaxation dynamics 
towards thermal equilibrium. While corresponding rates 
are of great interest experimentally, to extract them from 
PIMC data, typically requires very expensive long times 
simulations. 

Here we focus on this latter regime, not accessible by 
master equations, and develop a PIMC scheme which 
exploits the fact that a semi-group property applies ap- 
proximately. The proposed scheme works for P z {t) as 
observable as follows. One defines auxiliary densities 
p n {t) — Trs{W„(£)} on segments t 6 [nt s , (n + l)t s ] ,n = 
0, 1,2, ...N, where 

W n (t) = U{t,nt s ) Pn (nt s )^e-^ HB -^ U\t,nt s ). 

(18) 

and t s is sufficiently large compared to the bath memory 
time. In 

Pn(nt s ) = ~(l+ Yj p ^n(nt s )aA ,n>l, (19) 

\ v=x,y,z / 

the values for P u ,n(nt s ) as well as the bath-shift p n are 
adjusted such that one achieves a smooth matching ac- 
cording to 

Pz,n(nt s ) = P z , n -i{nt s ) , P z ,n(nt s ) = P e ,n-i(nt a ) , 
P z .A nt s) = P z .n-i{nt s ). (20) 

For n — 0, one has the initial density po(0) = | + 
with p,Q = 1. Now, due to continuity and the ex- 
act relation P y (t) = —P z (t)/A, the first two condi- 
tions fix the elements P z , n {nt s ) and P y n (nt s ) through 
the time evolution over the previous segment. To sat- 
isfy the third condition one properly adjusts P x ,n{nt a ) = 
5P x ,n + P Xi n-i(nt s ) and the bath shift p n . The reservoir 
distribution is thus progressively distorted by the spin 
dynamics and the reduced density carries information 
about how the bath modes induce the spin relaxation. 
This way, the correct dynamics over the total time inter- 
val ttot = {N + l)t s is approximated by a chain of PIMC 
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t=0.0, a=0.60, t s =1 5 
T=0.5, a=0.60, t s =1 5 
T=1.0, a=0.60, t s =8.0 
T=2.0, a=0.40, t.=8.0 




FIG. 2: Dynamics of P z (t) according to the CH-PIMC with 
optimized matching for various coupling strengths and tem- 
peratures. For each parameter set a different time interval t 3 
was used. See text for details. Temperatures are scaled in 
units of hA/kB- 



simulations over time intervals t s through 
P z (t) ~ P z , n (t) forte [nt a ,(n + l)t a ], n = Q,l,...N. 

Typical lengths of the segments range from a few up to 
15 times 1/A. The validity of this Chain-PIMC (CH- 
PIMC) is well-controlled. If for the combined data one 
obtains a smooth matching it works, if "kinks" appear, it 
does not. In the moderate to strong friction regime, this 
implies P x <c 1 which can be verified at each intermediate 
time step (for details see below) . Note that the procedure 
is exact at a — and thus expected to cover also the 
domain of very weak coupling. The computational gain 
is substantial: As the dynamical sign problem increases 
exponentially with t, the CPU time for a simulation over 
the total time interval t = nt s scales as e ants (a > is a 
constant) while for a chain with n segments of length t s 
one has an algebraic increase n e ats . 

In Fig.[5]results for the CH-PIMC are shown for Ohmic 
friction and larger couplings and over a wide tempera- 
ture range. Here, we found within the statistical errors 
an optimal matching for fi n — P z ^ n (nt s ) and SP x . n ~ 0. 
Hence, at each intermediate time step the drag of the sys- 
tem onto the equilibrium position of the bath is updated 
which accounts for dynamically induced system-bath cor- 
relations. The procedure also works in the domain of 
very low temperatures where the relaxation dynamics is 
strongly influenced by quantum effects but only as long as 
P x (nt s ) < 0.1. At zero temperature off-diagonal elements 
tend to be considerably larger, namely, P x {nt s ) « 0.3, 
and "kinks" at the intermediate time steps cannot be re- 
moved by varying 8P Xt n and In the domain Ah/3 < 2 
and a > 0.5 (lower temperatures, moderate to strong 
damping) one has thus access to the full equilibration 
dynamics where approximate treatments fail such as e.g. 
weak coupling master equations. When the reduced dy- 
namics is determined by strong coherences, but friction is 
not extremely weak, the CH-PIMC fails as seen in Fig. [3] 
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FIG. 3: Dynamics for P z (t) in the coherent regime (a — 0.3) 
using the CH-PIMC for various deviations 8P X , 8P y . Param- 
eters are T = and t 3 = 6/A. 



We recall that the CH-PIMC describes the bare coherent 
dynamics (a = 0) exactly though. 

VI. SUB-OHMIC HEAT BATH 

The methodology developed in the previous section 
is now used to analyze the SBM with sub-Ohmic spec- 
tral densities, i.e., densities ([FT]) with spectral exponents 
< s < 1. This class of reservoirs has received con- 
siderable attention recently as it seems to constitute the 
dominant noise source in solid state devices at low tem- 
peratures such as superconducting qubitsiS and quan- 
tum dots2£ with the spectral exponent s determined by 
the microscopic nature of environmental degrees of free- 
dom. It also appears in the context of quantum impurity 
systems^ and nanomechanical oscillators^. 

The model is further of fundamental interest with 
respect to its thermodynamical as well as its non- 
equilibrium properties. In thermal equilibrium at zero 
temperature a dissipation induced quantum phase tran- 
sition occurs^—. Accordingly, at critical coupling 
strengths a c (s) the system exhibits a transition from a 
delocalized state [P z (t —> oo) = 0] with tunneling be- 
tween the two spin orientations (weak friction) to a lo- 
calized one [P z (t — > oo) ^ 0] with basically classical be- 
havior (stronger friction) occurs. With respect to the 
relaxation dynamics, it has been shown recently^ that 
there is a coherent-incoherent transition in the regime 
1/2 < s < 1 at coupling strengths aci(s) < a c (s), while 
such a transition is absent with coherences surviving even 
for ultra-strong coupling in the domain < s < 1 /2 
(cf. Fig. In this regime, damped oscillations occur in 
P z (t) with frequencies that are in agreement with those 
extracted from the NIBA 



fi, w A., 



(22) 



Here A s = / °° du> J s {uS) / '(nu>) = 2aw c |r(s — 1)| is the re- 
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FIG. 4: Domains of coherent (white) and incoherent (shaded) 
dynamics of the SB model for a sub-Ohmic environment (|lip 
at zero temperature. Above (below) the solid line a c (s) the 
system asymptotically reaches a thermal equilibrium which is 
localized (delocalized)^. A coherent-incoherent changeover 
only occurs along the dashed line aci(s)^. 



organization energy. Likewise, the damping rate is given 
by 7 S (a) w sA s such that i} s /j s 1 (underdamped 
regime) . 

In this section we will explore first to what extent co- 
herences are robust against finite temperatures and var- 
ious cut-off procedures to better understand the domain 
in parameter space where the CH-PIMC is expected to 
work. Second, in the so determined incoherent regime 
relaxation rates are obtained from CH-PIMC data and 
compared to NIBA predictions. 



A. Coherent regime 

The central ingredient in the influence functional (j8]) 
which captures the impact of the reservoir is the bath 
correlation ([§]) . For sub-Ohmic spectral densities of the 
form (fTT|) it can be calculated analytically to read Q = 
Qo + Qp with the zero temperature part Qo(t) — 2aT(s — 
1)[1 — where z = 1 + ioj c t, and a part for finite 

temperatures 

Q fj {t) = 2ar(s-l> 1 - s [2C(s-l,l + l/K) 

- C(s -1,1 + z/k) - C(s -1,1 + z*/k)] (23) 

with the Hurwitz function C(z,a) and k = Tif3oj c . 

Now, according to the NIBA expression (fT2")l with Q 1 = 
Q'o + Qp an d Q" — Qq, finite temperatures lead to an 
exponential suppression of the integral kernel and have 
thus the tendency to destroy coherences (oscillations) in 
the relaxation dynamics of P z (t)- In particular, in the 
domain s <C 1 and for temperatures uj c Til3 ^> 1, one has 
Q'(t) ss Q' a {t)[l + t/(shf3)] so that on the time scale of 
the bare dynamics 1/ A thermal fluctuations are relevant 
if Ah(3 < 1/s. Then, due to the accumulation of low 
frequency bath modes, the SBM becomes progressively 
sensitive to thermal fluctuations for smaller s. 




FIG. 5: Dynamics of P z for various temperatures for s = 1/4 
and q = 0.1. The initial bath preparation is -Pz(O) = 1. 



This is illustrated in Fig. [5] which displays the 
changeover from oscillatory to overdamped motion at 
s = 1/4 with increasing temperature. Indeed, for tem- 
peratures fceT ~ shA coherences are washed out and 
after a transient period of time the dynamics turns into 
a monotone relaxation described by a single time con- 
stant, the relaxation rate. It is this rate that we extract 
below within the CH-PIMC. 

Before we do so, however, let us analyze the sensitiv- 
ity of the coherent dynamics on the cut-off procedure. 
According to (I22[) it is mainly determined by the reorga- 
nization energy A s . The corresponding integral is domi- 
nated by the low frequency profile of the spectral density 
J s {oj). As a consequence, for sufficiently large cut-off fre- 
quencies w c /A 3> 1 dynamical properties of P z (t) are 
basically independent of the specific form of the cut-off, 
e.g., a smooth exponential as in (|11[) or a sharp cut-off 
as in J s (uj) = 2iraujl~ s u! s 0(lo c — ui). This is in contrast 
to an infrared cut-off u>i which has pronounced effect on 
coherences as seen in Fig. [6] for T = 0. Even a small 
loi/A — 0.05 pushes the dynamics toward an incoher- 
ent decay similar to that at elevated temperatures, while 
varying the ultra-violet cut-off scheme has almost no in- 
fluence but only leads to a renormalization of the tunnel 
coupling 1 . To summarize this analysis, we find incoher- 
ent relaxation dynamics in the sub-Ohmic SBM for finite 
temperatures also in the domain < s < 1/2. The same 
is true for a low frequency cut-off, while the dynamics is 
robust against various high cut-off schemes. 



B. Incoherent regime 

In this regime, we start as above for the Ohmic SBM 
with the dynamics of P x (t) in Fig. [7] for various tem- 
peratures and coupling strengths. As expected, stronger 
friction and finite temperatures lead to a suppression of 
off-diagonal elements of the reduced density and thus 
drive the system towards a classical distribution. The 
NIBA expression (jl"5j) captures the PIMC data qualita- 
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FIG. 6: Dynamics of P z (i) for various cut-off schemes: a sharp 
cut-off (red, green, blue) and a soft cut-off (light blue, yellow, 
black) for the ultra-violet cut-off are shown together with the 
situation of an additional infrared cut-off (magenta); see text 
for details. 




FIG. 8: Polarization P z (t) for sub-Ohmic spectral density at 
T = hA/k B with s = 3/4, a = 0.5 and s = 1/4, a = 0.1 (see 
also Fig. [4]) using the CH-PIMC for various deviations SP X , 
5P V . 



T=0, 01=0.16 
T=0.1 , <x=0.20 
T=0.3, a=0.30 
T=1 .0, a=0.50 
T=2.0, a=1.00 




FIG. 7: Dynamics of the coherences P x (t) in a sub-Ohmic 
bath with s — 1/2 for various temperatures and coupling 
strengths. Symbols refer to PIMC data whereas solid lines 
show the corresponding NIB A results according to (| 15f) . 



tively correctly, but quantitative agreement is only found 
for sufficiently high temperatures and stronger couplings. 
In the low temperatures range, the NIBA predictions pro- 
vide a reliable estimate only for long times. Further, ac- 
cording to Fig. S] and the data in Fig. [7J one expects the 
CH-PIMC to be applicable for s > 0.5 and sufficiently 
high temperature, i.e. Ah/3 < 2. 

This is illustrated in Fig. [51 where the relaxation dy- 
namics for the polarization P z (t) is shown. Due to the fi- 
nite temperature, P z (t) decays in both cases monotonely 
but in the regime of small s- values, where at T = coher- 
ences survive up to arbitrary strong coupling (cf. Fig. 2]), 
off-diagonal elements are still substantial at weak fric- 
tion. One also finds that for a sub-Ohmic bath optimal 
matching is achieved for fj, n = 1, i.e. the initial bath 
preparation. The reason for that is the sluggish bath 
with a large portion of low frequency bath modes that 
remain basically static on the time scale of a few 1/A. 



In the range of spectral densities where the CH-PIMC 
applies, it can now be used to extract from the long time 
dynamics relaxation rates. However, as the required time 
disctretization q/t is much smaller for low temperatures 
and weak couplings, standard PIMC results are already 
sufficient to fit the rates in the regime where CH-PIMC 
usually fails, i.e. Ah/3 > 2. Our results can also be 
used to test the accuracy of the NIBA [cf. (|14l) ]. In fact, 
an analytical evaluation 14 yields the compact expression 
Tniba = «s exp [— 2ab s T(s — 1)] with 



b s = l + 2( Wc ?i/3) 1 - s [C( S -l,l)-C(s- 1,1/2)] 



2Ur 



ir(uj c hl3) 



l + s 



2oT(l + s)C(l + *, 1/2) 



1/2 



(24) 



Apart from the rather weak a-dependence in the pre- 
factor, the decay rate depends exponentially on the cou- 
pling strength. This is depicted in Fig. [S]with an almost 
linear behavior of log(l/r) as a function of a for dif- 
ferent temperatures. The NIBA gives indeed a perfect 
description over a broad range of parameters. The over- 
all magnitude of the lifetimes and thus the time scale for 
equilibration is extremely large compared to the bare dy- 
namics (see alsoii). Effectively, the spin is then almost 
frozen and displays a slow equilibration which may be 
beyond typical experimental time scales. These results 
reveal the benefit of the CH-PIMC, namely, to capture 
the domain of long times and strong coupling. One may 
argue that we test the validity of an approximate analyti- 
cal treatment (NIBA) by comparing its predications with 
those from an approximate numerical treatment (CH- 
PIMC). To avoid any inconsistencies we thus also per- 
formed full PIMC simulations and compared them with 
CH-PIMC results (not shown). The latter are again in 
perfect agreement within the statistical error with exact 
data in all cases where the matching procedure yields 
smooth relaxation dynamics. 
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FIG. 9: Inverse equilibration rates (lifetimes) in the incoher- 
ent regime for s — 1/2 as a function of a and T (in units 
of hA/ks)- NIB A predictions (lines) according to (|24[) are 
shown together with data extracted from standard PIMC 
(T < 0.3) and CH-PIMC (T > 0.5) simulations. 



C. Entropy and phase transition 

With the full reduced density operator at hand, an 
interesting quantity to measure with PIMC simulations 
is the von Neumann entropy 

S/k-B = -Tr(plogp) = -p+ \ogp+ - p- \ogp- (25) 

where p± = \ (l ± ^{a x ) 2 + (a y ) 2 + (cr z ) 2 ) . It is 

known to provide at zero temperature the entanglement 
between system and reservoir for the SBM2 4 -. In partic- 
ular, in case of the sub-Ohmic SBM it was shown pre- 
viously that at T = the quantum phase transition in 
thermal equilibrium from the delocalized to the local- 
ized phase corresponds to a maximum of the entrop y 24 ' 27 . 
Here we use PIMC techniques to explore to what extent 
dynamical information on finite time scales allows to de- 
tect this transition. 

We start in Fig. [TU] with a spectral density with ex- 
ponent s = 3/4. Apparently, towards the end of the 
simulation time range the entropy tends to saturate for 
increasing coupling strength. The PIMC simulations 
are not fully equilibrated, but from the data one can 
at least estimate a maximum in the entropy to occur 
for 0.3 > a > 0.24 with the correct valued being 
a c (s — 0.75) = 0.30. For longer simulation times the 
dynamical sign problem drastically deteriorates the data. 

This problem becomes more severe for smaller expo- 
nents s < 1/2 where PIMC simulations become increas- 
ingly demanding due to longer saturation time scales. 
In order to reduce the transient regime, we start not 
with a polarized initial state but from a density po(0) 
with elements P„ chosen closer to known thermal equi- 
librium values^, i.e., \i = P z (0) = 0.6, P x (0) = 0.8 and 
P y (0) = 0. As seen for s = 0.25 in Fig. [TTJ this proce- 
dure leads on shorter time scales to a maximum of the 
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FIG. 10: Entropy for various couplings a at s = 0.75 and 
T — 0. The inset shows a blow-up for longer times. See text 
for details. 

entropy within the range 0.02 < a < 0.025. This is in full 
agreement with the exact value a c (s = 0.25) rs 0.022. 




FIG. 11: Same as in Fig. [I0]but for s = 0.25 with an initial 
preparation P x (0) = 0.8, fx = P*(0) = 0.6, P„(0) = 0. 



VII. CONCLUSIONS 

In this work the time evolution of the reduced density 
of a spin-boson model has been investigated within nu- 
merically exact PIMC simulations in domains of strong 
coupling to environmental degrees of freedom and also at 
elevated temperatures. An efficient formulation has been 
developed to extract simultaneously information about 
populations and coherences within a single PIMC run. 
With the full density matrix at hand, segments of PIMC 
simulations on moderate time scales could be combined 
based on successively determined initial states. This 
method (CH-PIMC) allows to cover long times and to ex- 
tract e.g. equilibration rates which are not accessible by 
perturbative treatments. Specific applications for spin- 
boson models with sub-Ohmic spectral densities reveal 
the domains in parameter space where the method is re- 
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liable. Where it does not work, it gives insight into the 
strength of long-ranged system-bath correlations. 
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Appendix: Implementation of the operator for 
coherences 



The PIMC is only efficient when one extracts observ- 
ables at intermediate times < u from a single run over 
a total time t > u. Since the influence functional is non- 
local in time, the measurement operator in the PIMC 
corresponding to a x includes also contributions from the 
bath. One starts with a Trotter-Suzuki discretization of 
the time axis with q Trotter steps according to t = qr 
with time increment r. The discretized form of the influ- 
ence functional (jHJ in the original forward and backward 
paths a and a' reads^S 



1 3 

* = EoE 

3=2 k=2 


q q j— l 

A^fcOb + ci E A ? } ti + E o E 

3=2 j=2 k=2 


^^Af^-ME^O (A.l) 

7=2 7=1 


with 




A = 


Re [Q(t)], 


Al<i<g-2 
Ai<j< g _2 


Re [Q((j + l)r) - 2Q(jr) 
Im +Q((j - l)r)] 




Re [Q(r/2)], 


A (1) 

A 2<j<g _ 
A 2<3<q 


Re [Q((j - l/2)r) - Q((j - l)r) 
Im -Q((j-3/2)r)+Q((i-2)r)], 


= 


-Adm [Q(r/2)], 


^2<3<q - 


-Mlm [Q((j - l/2)r) - Q((j - 3/2)r)]. 



(A.2) 

Since paths are closed at the end of the full simulation 
period t one has Cg+i = 0. Further, according to the 
chosen path discretization, the spin flips occur at tj = 

(i-D- 

Now, let us first explain the procedure for vanishing 
coupling to the bath a — 0. Populations P\ = P z ((i — 
1)t) are then simply obtained from 



PI = 



E 



n k * ' a j+i) K ( a 3+i> °3 

3 = 1 



E 

i+l 



11 K * ( a 3 ' a 3 + l) K ( a 3 + l > 



x<5, 



(A.4) 



Accordingly, at each intermediate time step can be 
read out from a Monte Carlo run over the whole time 
interval < t' < t by applying the operator M"~ a = 

For off-diagonal elements P x = P x ((i — 1)t) things be- 
come a bit more complicated. To extend the expression 
for intermediate times 



(-1 



n = E 

<Tl,...,c 

XP0(c r l,O'l)^<7 < ,o-<(^o-i,l - ^O-i.-l) 



11 A '* (°j ' ^+1)^(^ + 1 » a 3 ) 



(A.5) 



such that one can extract expectation values of <r z and a x 
from a single PIMC run, one has to propagate artificially 
to a diagonal state at i + l. Thus, one multiplies (|A.5|) 
with 



1 9 E "Wi.^i+i 



CT i + l> CT i 



K*(p' i ,o J iJr i)K(p i + u a i ) 



(A.6) 



so that the measurement operator reads 



M a=0 = 



2A-*(^,< fl )Jf( - mj0 - i 



(A.7) 



(A.3) 



In the dissipative case a^Oan additional difficulty arises 
due to the specific path discretization which does not oc- 
cur for P 2 -observables. Namely, the naive ansatz to use 
M"j in the expression containing the influence func- 
tional is not correct. To measure P*, one also has to 
include interactions of Q (difference paths) with spins 
at earlier times j < i due to the non-locality in time 
of the influence functional. It is important to realize 
that due to the discretization the value of Q corresponds 
to the value of the spin path in the time interval 
(i — |) t < t < (a - |) t so that one has to drop the 
contribution of the influence functional stemming from 
the second half of this interval (i — 1)t . . . fi — i) t. To 
extract expectation values for a x and a z from the same 
PIMC sampling, the operator for off-diagonal elements 
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must thus include a 'compensation factor' exp R&° ff ] with where 



r(i-i)r 

= / dt'QO (t' - (i - 1) 

Jo 



$f - / di'Ci6 (*' - (i - l)r) 

r*' 



dt" - t")C{t") + iL"(t' - t")r](t")] 

,(i-*)r Af - Re [Q(r) - Q(r/2)] 

- i/x / dt'QQ (f - (i - l)r) Q"(i') . (A.8) A?!^,.! Re [Q((j - l/2)r) + Q((j + l)r) 

° Af Im -Q((j + l/2)r) + Q(jr)] 

Here, m is due to an initially shifted bath prepara- A off,i Re [Q(( . _ 1/2)t) _ 2Q(( . _ 1)r) 

tion. ihis correction factor vanishes for diagonal ele- ^off.i = j m +Q((i — 3 /2)r)l 

ments (Ci = 0). This way, one has for the measure- « 

ment of P*(u),0 < u < t in the PIMC the operator jf(f).off = _ vlm _ i / 2 ) r ) - - 3/2)r)] . 

M Xji = M°7° exp($° ff ) with the discretized form ' ^ A 

i i—1 

i=i i=i 



i •> -i^off.l . j- «off.l/- . -{>(^i),off 
+1&X: 771 + Q A » Cl + l(i X i : 



(A.9) 
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